/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"

#include "uhm/matrix/uhm/matrix.hxx"
#include "uhm/matrix/uhm/fla.hxx"

namespace uhm {

  static inline int solve_chol_1_flat( int fs, int ss,
				       linal::Flat_ ATL, linal::Flat_ ABL, 
				       linal::Flat_ t,  linal::Flat_ b );
  
  static inline int solve_chol_1_hier( int fs, int ss,
				       linal::Hier_ ATL, linal::Hier_ ABL, 
				       linal::Hier_ t,  linal::Hier_ b );
  
  static inline int solve_chol_2_flat( int fs, int ss,
				       linal::Flat_ ATL, linal::Flat_ ABL, 
				       linal::Flat_ t,  linal::Flat_ b );
  
  static inline int solve_chol_2_hier( int fs, int ss,
				       linal::Hier_ ATL, linal::Hier_ ABL, 
				       linal::Hier_ t,  linal::Hier_ b );
  
  
  void Matrix_FLA_::solve_chol_1_x() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix 
    // ----------------------------------------------------------
    solve_chol_1_hier( this->fs, this->ss,
		       this->hier.ATL, this->hier.ABL,
		       this->hier.xt,  this->hier.xb );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_chol_1_flat( this->fs, this->ss,
		       this->flat.ATL, this->flat.ABL,
		       this->flat.xt,  this->flat.xb );
    
#endif
  }

  void Matrix_FLA_::solve_chol_2_x() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    solve_chol_2_hier( this->fs, this->ss,
		       this->hier.ATL, this->hier.ABL,
		       this->hier.xt,  this->hier.xb);
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_chol_2_flat( this->fs, this->ss,
		       this->flat.ATL, this->flat.ABL,
		       this->flat.xt,  this->flat.xb);
#endif
  }    
  
  void Matrix_FLA_::solve_chol_1_r() {
#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix 
    // ----------------------------------------------------------
    solve_chol_1_hier( this->fs, this->ss,
		       this->hier.ATL, this->hier.ABL,
		       this->hier.rt,  this->hier.rb );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_chol_1_flat( this->fs, this->ss,
		       this->flat.ATL, this->flat.ABL,
		       this->flat.rt,  this->flat.rb );
    
#endif
  }

  void Matrix_FLA_::solve_chol_2_r() {
#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    solve_chol_2_hier( this->fs, this->ss,
		       this->hier.ATL, this->hier.ABL,
		       this->hier.rt,  this->hier.rb);
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_chol_2_flat( this->fs, this->ss,
		       this->flat.ATL, this->flat.ABL,
		       this->flat.rt,  this->flat.rb);
#endif
  }    

  static inline int solve_chol_1_flat( int fs, int ss,
				       linal::Flat_ ATL, linal::Flat_ ABL, 
				       linal::Flat_ t,  linal::Flat_ b ) {
    if (fs)
      FLA_Trsm( FLA_LEFT, FLA_LOWER_TRIANGULAR,
		FLA_NO_TRANSPOSE,
		FLA_NONUNIT_DIAG, FLA_ONE,
		~ATL, ~t );

    if (fs && ss)
      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
		FLA_MINUS_ONE, ~ABL, ~t,
		FLA_ONE, ~b );

    return true;
  }

  static inline int solve_chol_1_hier( int fs, int ss,
				       linal::Hier_ ATL, linal::Hier_ ABL, 
				       linal::Hier_ t,  linal::Hier_ b ) {
    
    if (fs)
      linal::trsm( FLA_LEFT, FLA_LOWER_TRIANGULAR,
                        FLA_NO_TRANSPOSE,
                        FLA_NONUNIT_DIAG, FLA_ONE,
                        ATL, t );

    if (fs && ss) 
      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                        FLA_MINUS_ONE, ABL, t,
                        FLA_ONE, b);

    return true;
  }

  static inline int solve_chol_2_flat( int fs, int ss,
				       linal::Flat_ ATL, linal::Flat_ ABL, 
				       linal::Flat_ t,  linal::Flat_ b ) {
    if (ss && fs) 
      FLA_Gemm( FLA_CONJ_TRANSPOSE, FLA_NO_TRANSPOSE,
		FLA_MINUS_ONE, ~ABL, ~b,
		FLA_ONE, ~t );
  
    if (fs) 
      FLA_Trsm( FLA_LEFT, FLA_LOWER_TRIANGULAR,
		FLA_CONJ_TRANSPOSE,
		FLA_NONUNIT_DIAG, FLA_ONE,
		~ATL, ~t );

    return true;
  }

  static inline int solve_chol_2_hier( int fs, int ss,
				       linal::Hier_ ATL, linal::Hier_ ABL, 
				       linal::Hier_ t,  linal::Hier_ b ) {
    if (ss && fs) {
      linal::gemm( FLA_CONJ_TRANSPOSE, FLA_NO_TRANSPOSE,
                        FLA_MINUS_ONE, ABL, b,
                        FLA_ONE, t );
    }

    if (fs)
      linal::trsm( FLA_LEFT, FLA_LOWER_TRIANGULAR,
                        FLA_CONJ_TRANSPOSE,
                        FLA_NONUNIT_DIAG, FLA_ONE,
                        ATL, t );

    return true;
  }
}
